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ABSTRACT 



The optimization of spacecraft trajectories in vacuum has 
received extensive consideration since the inception of space 
flight, yet, the effects of atmosphere have been largely 
neglected. The advent of low Earth-orbiting, large satellites 
and platforms necessitates that atmosphere be included in the 
optimization process. A practical means of studying this topic 
is as a problem in minimum-fuel orbital maintenance. Optimal 
control theory advances the notion that orbital maintenance is 
optimized through periodic thrusting as opposed to forcing 
Keplarian motion by nullifying the effects of drag with 
thrust. Further, this must be optimized by primer vectoring. 
This thesis examined the efficiency of a simple method of 
orbital maintenance using fixed-angle transverse thrusting. 
Results show that for the purpose of fuel-minimization, the 
width of the radial band in which the satellite is to be 
maintained, is dependent upon thruster size. In nearly all 
cases, a thrust-angle of 70 degrees maximized the fuel saved. 
This thesis shows that fixed-angle transverse thrusting does 
not improve on forced Keplarian motion and hence thrust 
vectoring must be optimized. 
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NOMENCLATURE 



a r radial acceleration 

a tr transverse acceleration 

a semi-major axis of an ellipse 

C negative reciprocal of thruster exhaust velocity 

C D coefficient of drag 

D magnitude of drag 

D r radial component of drag 

D tr transverse component of drag 

E total energy 

F r external forces in radial direction 

F tr external forces in transverse direction 

g G gravitational acceleration at sea level 

H Pontryagin's H-function 

I sp specific impulse 

f_ unit vector in direction of thrust 

m spacecraft mass 

r_ position vector 

r radius 

r rax maximum radius 

r min minimum radius 

R 0 initial radius 

R f final radius 

S path travelled by spacecraft 

51 initial position of spacecraft 

52 final position of spacecraft 

S ref aerodynamic reference surface of vehicle 

T thrust vector 

T„. ax magnitude of thrust 

T r radial component of thrust 

T tr transverse component of thrust 

t time 

t Q initial time 



vii 



]>> 1 ^ 

N < 



t f final time 

v velocity 

v velocity vector 

v, radial component of velocity 

v t , transverse component of velocity 

a thrust angle 

(3 exponential density scale factor 

AE change in total energy 

AR radius at which periodic thrusting is commenced 

At change in time 

% specific energy 

0 orbital coordinate used to define spacecraft position 

X Lagrange multiplier 

X, mass costate 

velocity costate vector 
radius costate vector 
(1 gravitational constant 

p atmospheric density 

p 0 atmospheric reference density 
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I . INTRODUCTION 



In 1963, D.F. Lawden, in his monograph [Ref. 1], 
laid the foundation for what has become a highly sophisticated 
subdiscipline of astrodynamics , optimization of space 
trajectories. In 1979, Marec [Ref. 2] provided a more 
comprehensive treatment of the subject in his text on optimal 
space trajectories. Examination of the optimization of 
spacecraft trajectories has been treated by many authors in 
manners similar to these two great works, yet, until recently, 
the study of atmospheric effects upon these trajectories was 
largely neglected. Research and development of hypervelocity 
vehicles have kindled interest in this area, through which 
study, other areas of interest have emerged. 

One such area is the effect of aerodynamic force on non- 
lifting, or, blunt bodies. First addressed by Ross and Melton 
[Ref. 3], this subject is of particular interest for 
two reasons. First, atmospheric effects on low-Earth orbiting 
(LEO) satellites are of obvious interest and second, as stated 
in their paper [Ref. 3:p. 2], better understanding of this 
phenomenon could provide deeper insight into the more 
complicated topic of lifting bodies in the upper atmosphere. 

Much work has been accomplished in this area pertaining to 
accurate prediction of satellite orbits [Ref. 4, 5]. 
The focus on atmospheric effects as they pertain to the 
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specific problem of minimum-fuel orbital transfers, however, 
is unique to [Ref. 3] and the follow-on work described here. 
The particular problem of orbital maintenance can provide 
insight into the more general, and complex problem of orbital 
transfer . 

Historically, orbital transfers (coplanar) have been 
accomplished as either, two, or, three-impulse maneuvers 
[Ref. 6:pp 78-88]. The problem of orbital transfer is 
approached as a minimization of energy required to move a 
satellite from orbit A to orbit B, or, equivalently, a 
minimization of the characteristic velocity. A satellite that 
has descended from an initial orbit due to a disturbing force 
such as drag, and which must be returned to its initial orbit 
can be approached as a problem requiring an orbital transfer. 
Orbital transfers such as this are optimally performed by two- 
impulse transfers, such as a Hohman transfer. Large orbital 
transfers (r a >11.8r F ) are optimized with three-impulse 
maneuvers [Ref. 6:p 87]. As low-Earth orbits become more 
frequently utilized, deeper understanding of the effects of 
drag must be achieved in order that orbits, propulsion systems 
and costs are optimized. This is particularly applicable for 
large satellites, such as the proposed space station, that 
must remain in low orbits for extended periods of time. 

Additionally, more complex areas of study, such as that of 
lifting bodies in the upper atmosphere, could benefit from the 
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insights gained through a deeper understanding of atmospheric 
effects on non-lifting bodies. 

Research on the atmospheric effects on low-Earth orbiting 
spacecraft is sure to receive much attention in the future. 
The benefits to existing and future systems, while extensive, 
remain relatively unexplored and demand the attention of the 
aerospace industry. 
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GENERAL FORMULATION OF THE PROBLEM 



II . 

Optimization of orbital transfers is a subject that has 
achieved a high degree of sophistication and many elegant 
solutions exist [Ref. 1, 2]. However, the specific treatment 
of non-lifting bodies is in the initial stages of development. 
The increasing number of low Earth orbiting satellites 
requires that a study of atmospheric effects on orbital 
trajectories be conducted. In this thesis, the problem of 
minimum fuel orbital maintenance is considered. Two methods 
are examined by which orbital maintenance may be performed. 

One method is to counter drag with thrust. In this forced 

Keplarian motion, the reaction control system would thrust 

continuously for the duration of the satellite' s lifetime with 

magnitude and direction equal and opposite to drag. The second 

scheme considered here, utilizes periodic transverse 

thrusting, or, non-Keplarian motion to correct for 

perturbations due to drag. The question to be answered is 

whether or not optimal non-Keplarian trajectories are superior 

to f orced-Keplarian trajectories. 

/ 

Let the problem be defined as maintaining an orbital 
deviation within a specified radial band r mln <r<r max . Is forced- 
Keplarian orbital maintenance, i.e., exactly countering 
aerodynamic forces with thrust, superior to non-Keplarian 
orbital maintenance, i.e., allowing the orbit to decay to a 
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certain point, then reboosting to a point above the desired 
altitude. While maintaining the spacecraft within the 
specified radial band? The study performed by Ross and Melton 
[Ref. 3:p. 4] suggests that f orced-Keplarian motion is not 
optimal and that thrust vectoring must be considered if an 
optimal solution is to be obtained. 

This thesis addresses an additional question, can orbital 
maintenance be optimized if thrusters are fired at a fixed 
angle to the local horizon and if so, what is the angle or, 
preferably, range of angles that achieve optimality? 

Ross and Melton [Ref. 3] develop their theory through the 
methods of optimal control theory. It is proposed here that, 
while ideally accomplished in that manner, satisfactory and 
enlightening results may be obtained through the use of 
relatively unsophisticated mathematics and the aid of computer 
modelling techniques. To exemplify this statement, let us 
examine the mathematical development in Ross [Ref. 3: pp. 
1-3] . Drag is given by 

D = ~±P(Z) S zef C D vy ( 1 ) 

where p (r) is atmospheric density, S ref is the reference 
surface area of the spacecraft upon which the aerodynamic 
forces act, C D is spacecraft's coefficient of drag and v is 
its velocity. The equations of motion can be written as 

t=Y ( 2 ) 
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^zef^D 



( 3 ) 



v = 




fLXv + 

m m 



Applying the 
given by 



m = cr M ,; 



C = — 



SJ & o 



principles of Pontryagin, the 



( 4 ) 

Hamiltonian is 



H - A r v - 





T 

m — ^ 



A m cr m 



( 5 ) 



The costates are then written as 




i _ dH 
m dm 





( 8 ) 



A closed form solution to these equations does not exist and 
only after initial and boundary conditions have been 
determined may the solution be obtained through numerical 
integration. This is a cumbersome method for determining the 
optimal direction of thruster firing. In this thesis, we look 
into the possibility of a constant vectoring scheme that may 
result in nearly identical performance. 
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A computer generated model will strive to maintain 
spacecraft trajectory within ±1.5 kilometers of the injection 
altitude utilizing a periodic fixed angle transverse thrusting 
control scheme whose direction is maintained at a fixed angle 
relative to the local horizon. The propellant consumed will 
then be compared to that of the same satellite employing a 
control law that sets thrust equal to drag at every point 
within the orbit. 



7 



III. SIMPLIFIED PARAMETRIC FORMULATION 



A method less elegant than optimal control theory, but 
nonetheless valid, is that of parametric variation. The 
equations of motion are developed for the satellite's orbital 
trajectory and certain parameters varied to achieve "optimal” 
control of orbital variations. 

A. DEVELOPMENT OF THE EQUATIONS OF MOTION 

For the purpose of building a solid foundation, certain 
simplifying assumptions are made. Orbital motion is assumed to 
be coplanar, the initial spacecraft orbit is assumed to be 
circular and since the spacecraft is a non-lifting (blunt) 
body, drag is the net aerodynamic force acting upon it. 

The equations of motion for this two-body system can be 
written as 



a r 




( 9 ) 



a =£— (10) 

where a r is the radial acceleration of the spacecraft, a tr is 
its transverse acceleration, F r is the sum of the external 
forces in the radial direction, F tr is the sum of the external 
forces in the transverse direction and m is spacecraft mass. 
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The external forces acting on the spacecraft are the 
gravitational field, aerodynamic forces and thrust. Figure 1 



illustrates the coordinate system and the net forces. 




Figure 1 Graphical Representation of Coordinate System 

Referring to Figure 1, it is clearly seen that the components 
for drag are given by 



D r - -D sin (y) 


(11) 


D tz = -D cos (y) 


(12) 
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Likewise, thrust is written as 



T i = T n*x cos (a) 



(13) 



T tI = T t ^ sin (a) (14) 

The angle y is the flight path angle, defined by the 
intersection of the velocity vector and the transverse axis. 
The angle a is called the thrust angle, defined by the 
intersection of the thrust vector and the transverse axis. The 
equations of notion, then, become 



I - e 2 r = + II 


(15) 


z 2 iti m 




D T 

e r + 2 6r - - tr + tr 


(16) 



m m 



where |i is the Earth's gravitational constant. 



B. DETERMINATION OF THE CONTROL VARIABLE 

Following the development of the previous chapter, it is 
desirable to maintain the spacecraft within a radial band of 
a predetermined width. This may be accomplished by directly 
controlling either radius or eccentricity. Radius is the 
control variable of choice for a number of reasons: it is 
found directly from integration of the equations of motion, 
changes are easily visualized and radius control provides 
indirect control of the eccentricity. It is clear that by 
varying the thrust, control of satellite radius is possible. 
Examination of the thrust equations (Equations 13 and 14), 
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presents two methods by which thrust may vary, changes of 
amplitude or change in the direction of the thrust vector. 

The following chapter presents a method by which radial 
deviation is controlled through variation of the thrust angle 
and then tested by varying the thrust magnitude. A computer 
model is developed that simulates the trajectory of a 
spacecraft, graphically demonstrating the effects imposed on 
it through variation of the direction of the thrust. 



11 



IV. DEVELOPMENT AND TESTING OF THE COMPUTER MODEL 

A. COMPUTER PROGRAM DEVELOPMENT 

As indicated in the previous chapter, a computer program 
is developed to simulate spacecraft orbital trajectories and 
is included in Appendix A. The program is written in the 
fortran language and employs a fourth-order Runge-Kutta 
numerical integration routine to integrate the equations of 
motion. The program consists of six sections, the main program 
and five subroutines. The main program controls input and 
output while the subroutines provide various other functions 
necessary for accurate simulation of orbital trajectories. 

The first subroutine calculates drag experienced by the 
spacecraft. Initially, a model incorporating constant 
atmospheric density is used. This facilitates verification of 
the program, after which, an exponential density model is 
used. It is acknowledged that more accurate atmospheric 
density models exist, however, the exponential model provides 
satisfactory accuracy over the range of travel experienced by 
the satellite (± 1.5 kilometers from initial orbit R 0 ) as 
controlled by the simulation. The second subroutine contains 
the equations of motion governing the spacecraft's orbital 
trajectory. To facilitate handling, the equations were broken 
into parts. Solving Equations (15) and (16) for acceleration, 
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it is seen that the right-hand side of the radial equation has 
four components and the angular equation has three. These 
components are labeled A, B, C, and E for the radial equation 
and P, Q and S for the equation governing angular motion. The 
third subroutine contains the fourth-order Runge-Kutta 
numerical integration routine used to integrate the equations 
of motion. The next subroutine calculates the parameters of 
the satellite's osculating orbit. The last subroutine contains 
the control law governing the activation and deactivation of 
the thrusters responsible for the periodic maintenance of the 
satellite's orbit. 

The following parameters define the specifications around 
which the computer model was developed. 

• Spacecraft mass - 20,000 kg. 

• Spacecraft frontal area (S re5 ) = 60 m 2 . 

• Coefficient of drag (C D ) = 2.2. 

• Altitude of spacecraft's orbit (R c ) = 260 km. 

• Atmospheric density at R 0 (p) = 8.3130 x 10" 11 kg/m 3 . 

B . PROGRAM VALIDATION 

Program development proceeded in stages, with each stage 
requiring validation prior to beginning the next. Initially, 
all external forces except gravitation were neglected. 
Clearly, radius, speed, angular momentum and specific energy 
must remain constant for the program to be considered valid. 
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Having accomplished that elementary stage, phase two 
introduced drag (Equation (1)). Taking into account the fact 
that aerodynamic forces acting on a spacecraft are very small, 
and consequently, the time required for significant changes to 
occur, very large, atmospheric density was increased by 
approximately three orders of magnitude in order to reduce 
computer run time. 

Initially verifying that radius continually decreased, the 
accuracy was tested by comparing the difference in altitude 
per orbit calculated by the program to that calculated 
manually by simplified equations. This is accomplished through 
the use of a "rule-of-thumb" and is developed in the following 
section . 

1 . Development of a Rule of Thumb 

Work done by drag is a function of the path travelled 
by the spacecraft. Therefore, the amount of work done 
corresponds to the change in energy of the spacecraft, which 
is given by 

Work done = A E = Drag x 2nr (I 7 ) 
Spacecraft specific energy is given by 

r = zl - J± = -JL (18) 

2 r 2a 

Total energy E is equal to the specific energy multiplied by 
spacecraft mass. Assuming that the spacecraft is in a circular 
orbit, the semi-major axis a is equal to the radius r. 
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Performing this substitution then differentiating both sides 
of the equation while holding mass constant yields 

dE = -H£L dr (19) 

2r z 

Setting Equations (17) and (19) equal to each other and 
solving for the change in radius yields 

Ax = - 4 -?/ 3 - 0 (20) 

[irn 

Equation (20) is the decrease in radius per orbit of the 
spacecraft due to drag. Despite inaccuracies resulting from 
the simplifying assumptions, this rule of thumb is accurate to 
within a few percent. The difference in spacecraft radius 
calculated by the computer program matched that of the thumb 
rule within a few percent thereby validating the model through 
this point in its development. Phase three introduced thrust 
while setting drag equal to zero. Clearly, any results other 
than steadily increasing radius, angular momentum, and 
specific energy would have been cause for program 

invalidation . 

2. Development of a Control Law 

The purpose of the control law is to maintain 
satellite radius within a prespecified bandwidth. By 
monitoring certain variables, activation and deactivation of 
station keeping thrusters can be determined. Keeping the 
control law simple in order that computer memory and run time 
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related to station keeping be kept to a minimum is also a 
desirable goal. With these facts, control law design proceeded 
as described in the next section. 

3. Control Law design 

a. Approximate method 

Specific energy is a function of spacecraft radius. 
For that reason, specific energy is a very useful parameter 
that can be used to maintain the satellite within the 
specified bandwidth. The energy lost when the satellite's 
radius decreases, is dependent upon the path taken by the 
satellite in its descent as illustrated in Figure 2 below 
(exaggerated for clarity) . 




Figure 2 Path Travelled By A Spacecraft Between Two Orbits 

If S represents the path travelled by the spacecraft, then Sj 
is its initial position and S 2 its final position. The force 
acting against spacecraft motion is drag, which directly 
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opposes the velocity vector. With this in mind, the energy 
loss is given by 

A E = Jz? dS (21) 

Si 

The arc length S is a function of radius and angle turned 
through, and is given by 

5 = r6 (22) 

Differentiating the above equation yields 

dS = rdd + 6 dr (23) 

Substituting Equation (23) into Equation (21) yields energy 
loss in terms of the known variables, r and 0 

e 

A E = J Eftdr + jDrdd (24) 

*o o 

Employing Simpson's rule, the loss in energy can be 
approximated fairly closely. Assuming atmospheric density, 
flight path angle, thrust angle and thrust are all constant, 
Simpson's rule applied over ten iterations yields the energy 
lost in moving from R c to R ; 

ae = ("“ Sy "" 2 } 00 G0) + To Dr (29) (25) 

Referring to Figure 2, R f is the point at which the thrusters 
will fire. In terms of the control law, R f is the control 
variable and can be relabelled as AR since it is variable and 
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determines the width of band in which the satellite is 
maintained. As developed for this model, AR is the initial 
radius R 0 minus one kilometer. This value is chosen so as to 
maintain maximum orbital deviation within ±1.5 kilometers of 
R c . Now that the energy loss has been approximated, the 
objective is to determine the length of time to fire the 
thrusters in order to replace the energy. The return path of 
the spacecraft is a function of the thrust. 

A E = J (T - D) • ds (26) 

Recognizing that d_s is related to the time rate of change of 
the position vector or arc length _S, Equation (26) becomes 

A E = J (T - D)»v dT (27) 

Reducing Equation (27) to component form yields 



t[ c f 

A E = f ( T Z ~D Z ) v z dt + f (T cz - D tz ) v cz dt (28) 

Co Co 

Integrating and solving for At yields the length of time that 
the thrusters must fire in order to replace the energy lost 
due to drag. 



At = 



A E 



(29) 



(T z - D z ) v z + (T CI - D tz ) v tz 

Examination of Figure 1 reveals that velocity can be written 
as 
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v = vsin (y) 



(30) 



and 



v cr = vcos (y) (31) 

Figures 1 through 6 of Appendix B illustrate 
results obtained from this method for thrust angles of 60° and 
65°. Examination of plots of thrust versus orbit (Figures 2 
and 5) reveal the inadequacy of this model. Thruster firing 
times are seen to be on the order of orbits rather than 
fractions of orbits. This is due to two factors, first and 
most important being that the energy change is calculated from 
the initial orbit to the point AR where the thrusters begin 
firing. The problem arises from the fact that the satellite is 
still in a descent at this point and continues to descend 
until its motion is reversed through the opposing force of the 
thrusters. As a result, the satellite loses more energy than 
is replaced. The second problem arises from the inaccuracies 
inherent in the assumptions required to perform the 
approximation. While the first problem is easily resolved, the 
changes in atmospheric density and flight path angle, while 
very small, are not constant and the resulting errors combine 
to render this model unsatisfactory. A more accurate method is 
to calculate the energy loss directly using the variables 
derived from integrating the equations of motion. 
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b. Direct method 



Using Equation (18), the initial and instantaneous 
energies may be calculated at any instant during the 
trajectory of the spacecraft. As before, thrusters will fire 
when the spacecraft descends below AR. The program then 
calculates the spacecraft' s specific energy each time the 
equations of motion are integrated. Comparing this value to 
that of the initial circular orbit, the control law commands 
the thrusters to fire until they are equal. Examination of 
Figures 1 through 6 of Appendix C reveals that while results 
are closer to those expected, this method also appears to be 
inadequate. Thruster firing times (Figures 2 and 5) remain 
excessively long. Reevaluation of the model suggested that the 
solution might be a function of the thrust to drag ratio. 
Increasing thrust to a value of 300 N and then examining 
results for thrust angles of 60°, 65° and 70° revealed this to 
be the solution. Results are found in Appendix D. Thrust 
angles of 65° and 70° maintain radial deviation within the 
prescribed bandwidth of three kilometers with a trend that 
indicates they will remain so indefinitely. Thrust plots 
(Figures 2, 5 and 8) illustrate that the thrusters are firing 
over a small portion of an orbit instead over a period of many 
orbits as before. Plots of spacecraft radius versus orbit are 
included in each appendix to illustrate the success of the 
control law in maintaining radial deviation within the three 
kilometer band. Energy plots are also included as 



20 



corroberat ing illustrations of the spacecraft's energy level 
at each point in its trajectory. 
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V. ANALYSIS AND RESULTS 



Analysis of the validated model is performed in four 
steps. Prior to testing, drag was returned to its true value 
of 8.3130 x 10* :: kg/m 3 and thrust was set equal to a value of 
25 N. This value is determined by multiplying the thrust to 
drag value of the test model by the drag actually experienced 
by the satellite. The thrust to drag ratio with the increased 
density is 



3 00/7 
3.956 N 



75.834 



(32) 



With p = 8.313 x 10" :: kg/m 3 , drag is calculated to be 0.329 N. 
Multiplying this to the value in Equation (32) yields a thrust 
of 25 N. 



A. USE OF A CONSTANT DENSITY MODEL 

Initial testing maintained a constant density while 
varying thrust angle. The results obtained from these tests 
are interesting. As illustrated in Appendix E, approximately 
six and a half orbits are required for the spacecraft's radius 
to decay to AR. Thruster firing places the spacecraft into a 
slightly elliptical orbit (typically, eccentricities were 
found to be on the order of 10' 4 ) . Examination of energy plots 
reveal that the energy of the elliptical orbit is very close 
to that of the initial circular orbit. This results in one 
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large energy change as the spacecraft initially descends from 
R c and is returned, then, subsequent small changes are 
required once the spacecraft is established in its 
"elliptical" orbit. A thrust angle of 70° is seen to yield the 
best results (see Figure 6 of Appendix E) . Radial deviation 
remains just within the specified radial band with a trend 
that indicates it will do so indefinitely. At angles of less 
than 70°, radial deviation steadily increases until it exceeds 
the prescribed limits. An illustration is provided in Appendix 
E, for a thrust angle of 65°. Thrust angles greater than 70° 
follow a trend illustrated by Figure 7 of Appendix E (a=75°) 
until reaching approximately 85°. Above 85°, increasingly 
larger values of thrust are required to maintain the 
spacecraft within the prespecified radial band. 

B. USE OF AN EXPONENTIAL ATMOSPHERIC MODEL 

Use of a constant atmospheric density model permitted the 
determination of an optimum angle that could be compared with 
that determined by a more accurate exponential atmospheric 
density model. Also, realizing that density changes would be 
small within a three kilometer band, the constant density 
model divulged a reasonable facsimile of results obtained from 
the exponential model. As stated previously, more accurate 
density models exist, for example, the Jacchia Atmospheric 
Model (J70) . The exponential model, however, provides 



23 



sufficient accuracy for the development and testing done in 
this thesis. With this in mind, density is now given by 



P 



Po e 



-P <r-r 0 ) 



( 33 ) 



where |3 is determined to be -2.12 x 10 -5 m _1 [Ref. 7] . 
As predicted, the results are nearly identical to those 
obtained from the constant density model. A thrust angle of 
70° maintains the spacecraft within the three kilometer band 
with the same trends as described in the previous section. 
Results are illustrated in Appendix F. 

1. Comparison with a Forced-Keplarian Model 

It is now possible to compare results obtained from 
this model with that of a spacecraft experiencing forced- 
Keplarian motion (thrust equal to drag, resulting in a 
circular motion) . Modifications to the program to obtain a 
model in which thrust is equal to drag are very simple. The 
drag subroutine calculates drag then sets thrust equal to it. 
The subroutine containing the thrust control law is removed 
from the program entirely since thrusters will fire 
continuously as long as propellant is available. Drag always 
opposes the velocity vector which, as previously shown, is 
defined by the flight path angle y. If the thrust is exactly 
opposite and equal to the drag force, then the thrust angle a 
is equal to the flight path angle. Therefore, modifications 
consisted of setting a = y, removing the subroutine containing 
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the thrust control law and setting thrust equal to the drag 
calculated in the appropriate subroutine. Since the spacecraft 
is initially in a circular orbit, validation of these 
modifications is achieved when the spacecraft's orbital 
parameters remain unchanged over the test period; in this 
case, 20 orbits. 

Finding that the model performed as expected, plots of 
propellent consumed over the test period are compared with 
those for the spacecraft experiencing non-Keplarian motion at 
a thrust angle of 70°. Results are contained in Appendix G. It 
is clearly seen that orbital maintenance using forced- 
Keplarian motion is superior to that using non-Keplarian 
motion . 

2. Further Testing 

To determine the "robustness" of the model, two 
additional tests were performed. In the first, thrust was 
maintained at 25N while specific impulse was varied over the 
range of 200 sec, valid for chemical reaction engines, to 2000 
sec which is valid for electrothermal engines. In the second 
test, specific impulse is returned to its initial value of 300 
sec while thrust is varied over a range extending from IN to 
35N . 

a. Constant Thrust, Varied Specific Impulse 

Varying specific impulse while holding thrust 
constant, forces the rate of fuel consumption to change. 
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Specific impulse is given by given by 




where I sp is the specific impulse, T is the magnitude of the 
thrust, g 2 is the gravitational acceleration at sea-level and 
m is the change in mass. 

Clearly, decreasing thrust results in a corresponding 
decrease in the rate of propellant consumed over a given 
period of time. This is graphically represented in Appendix G. 
The absolute quantity of fuel consumed decreases with 
increasing specific impulse, but the percent difference 
between the reboost and forced Keplarian methods remain 
virtually unchanged.' 

The percent difference is calculated by taking the ratio 
of the values of propellant mass at a specific time for 
thrust-equal-drag plots and reboost plots to determine the 
relative change between them. This provides a truer comparison 
of the two schemes than does simply comparing the end values 
of the plots. 

The results illustrated in Appendix G indicate that, while 
playing a significant role in propulsion system optimization, 
specific impulse is not a function of the method used to 
perform orbital maintenance and will not affect the particular 
outcome of one method more than another. 
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b. Constant Specific Impulse, Varied Thrust 

The previous section illustrated the significance 
played by specific impulse in minimizing fuel expenditure 
during orbital maneuvers. Obviously, engine "size" plays an 
equally important role in that process. It is expected that as 
thruster size increases, fuel expenditure will increase. 
Appendix H bears this out. While it is theoretically a simple 
matter to choose the proper specific impulse (bigger is 
better), this is not the case when choosing thruster size. 
Examination of the plots in Appendix H reveal that when 
comparing engines over a certain time span, the end results do 
not provide a ready solution. Figure 1 is a case in point. 
Although this case (IN thruster, I sp =300 sec) results in the 
least amount of fuel expended over the range of thrusts 
chosen, it is obvious that this is not a wise choice of 
thruster size. The thruster fires continuously from its 
initiation until the end of the test period. The thruster is 
clearly too small. A five newton thruster, i.e.. Figure 2, 
seems to be a viable engine size, although, without more 
information, it is difficult to determine positively. 

Relatively large thrusters burn for shorter periods 
of time than smaller thrusters to achieve a common result, but 
each burn expends more fuel. Smaller thrusters, on the other 
hand, expend less fuel for a given burn time, yet must burn 
longer to achieve the same results. The fundamental result of 
this test is that thruster size is a tradeoff variable that is 
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to be used in conjunction with other factors to achieve 
desired results, such as minimization of thruster burn time 
during orbital maintenance. 

3. Final Tests 

Previous sections have shown that within a narrow 
radial bandwidth, a model using thrust equal to drag is 
superior to one using fixed thrust-angle reboost techniques. 
The question arises as to whether these results will remain 
valid for larger bandwidths . 

As described earlier, the control law commands the 
thrusters to fire when the spacecraft orbit has decayed a 
certain distance below the reference orbit. Results are 
examined for cases where the spacecraft is allowed to descend 
20 km, 30 km and 40 km below R 0 . Thrust is fixed at 300 N 
while specific impulses vary between 300 sec, 500 sec and 2000 
sec. As for the case of a 3 km bandwidth, a thrust angle of 
70° maintained the spacecraft within the desired bandwidth and 
was therefore used for all cases described below. The tests 
are performed over a period of 100 orbits. In order to reduce 
computer time, atmospheric density was once again increased to 
a value of 1 x 10' 9 kg/m 3 . 

a. Case 1: Thrusters Fire 20 km Below R 0 

Rather than choosing a specific bandwidth and 
adjusting the control law to achieve it, the spacecraft's 
orbit was allowed to decay a certain distance prior to 
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activation of the control thrusters and the resulting 
bandwidth measured. This provided expediency since the control 
law determined the bandwidth rather than having to be adjusted 
to achieve it. The results are the same in either case so no 
accuracy is lost with this method. The figures in Appendix I 
illustrate the results of this case. 

Allowing the spacecraft's orbit to decay 20 km 
prior to activation of the control jets provided a radial band 
of 55 km. Plots of expended propellant mass versus orbit are 
provided for the three test cases (I sp = 300 sec, 500 sec and 
2000 sec) . As in previous cases, thrust equal to drag provides 
a straight line while reboost is actually a series of steps. 
Each step is a cycle wherein the spacecraft descends 20 km at 
which time reboost occurs (vertical portion of plot) , after 
deactivation of the thrusters, the spacecraft once again 
descends to the point where reboost reoccurs. This is 
indicated by the horizontal portion of the plot since no fuel 
is being expended during this portion of the trajectory. As 
before, the reboost maneuver puts the spacecraft into a 
slightly eccentric orbit (on the order of lO -4 ) . This accounts 
for the periodic motion and high number of thruster firings 
indicated by the mass plots in Figures 2 through 4. As 
expected, increasing specific impulse reduces the amount of 
propellant expended over the test cycle. 
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b. Case 2: Thrusters Fire 30 km Below R 0 

Requiring the control law to fire 30 km below the 
reference orbit provides a radial band of 78 km. Figure 5 is 
an illustration of spacecraft radius over time. Because AR is 
so large, the reboost maneuver must occur twice before the 
spacecraft settles into a periodic trajectory that carries it 
the full width of the radial band. Eccentricity of the 
osculating orbit, however, remains on the order of 10' 4 . Once 
in its periodic trajectory, results are very similar to those 
of case 1. Figures 6 through 8 illustrate propellant mass 
consumption over the test period. The first two thruster 
firings are clearly evident. Once in its periodic trajectory, 
the mass plots are very similar to those of case 1 and occur 
for the same reasons. 

Comparisons of the mass plots of case 2 to case 1 
reveal interesting results. Percent difference comparisons of 
case 2 to case 1, for identical specific impulses, provides an 
indication of established trends from which inferences of 
future results might be made. 

We see that for case 1, for each variation of 
specific impulse, reboosting the spacecraft requires 482 
percent more fuel than the use of thrust equal to drag 
techniques. Case 2 reboost maneuvers require 578 percent more 
fuel than thrust-cancel-drag maneuvers. 
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c. Case 3: Thrusters Fire 40 km Below R 0 

Allowing the spacecraft's orbit to decay 40 km 
prior to initiating reboost sequences provides a radial band 
of 100 km in which the satellite trajectory is maintained. 
Figure 9 illustrates radial deviation of the spacecraft. The 
increase in AR coupled with a fixed thrust requires the 
satellite to perform three reboost maneuvers before the 
familiar periodic trajectory is achieved. Thruster firing is 
clearly evident in the first two incidences, as is the ensuing 
orbital decay of the resulting (slightly) eccentric 
trajectories . 

Calculating the mass expenditure percentages 
reveals that reboosting the satellite requires 481 percent 
more fuel than does setting thrust equal to drag. Although 
thrust equal to drag is still superior, the difference between 
the two is decreasing. To further test this result, 
percentages were calculated for points 43 orbits and 97 orbits 
into the test period. All values were less than corresponding 
values calculated in case 2. While these results do not 
provide conclusive evidence, we may conjecture that a trend is 
developing, indicating that at some point results from 
periodic orbital maintenance will equal those from forced 
Keplarian motion, or as the theory predicts, the plots will 
reverse and periodic thrusting will become superior. 
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C. SUMMARY OF RESULTS 



Initial testing, utilizing a constant atmospheric density 
model provided a baseline against which, further testing could 
be compared. A thrust angle of 70° was found to produce the 
desired results. Fixing the thrust vector at this angle 
maintained spacecraft orbital deviation within a three 
kilometer band nearly indefinitely. Upon determining this 
angle, the computer program discarded the constant density 
model and incorporated an exponential atmospheric density 
model . 

This simulation was then compared to a model in which 
thrust canceled drag. According to the optimal control theory 
developed in Chapter III, the simulation (reboost model) 
should prove superior to a thrust-cancel-drag model (relative 
to the problem of fuel-minimization) . In fact, reboost 
required significantly more propellant to maintain the 
satellite orbit within the three kilometer band than did 
thrust equal to drag. To test the robustness of these result, 
specific impulse was varied between 200 and 2000 sec and 
thrust was varied between 1 and' 35 N. The results remained 
unchanged. To further test the results, the simulation was 
tested over much wider radial bands. The results still proved 
thrust-cancel-drag trajectories superior to reboost models 
although the difference in efficiency seemed to decrease. 
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VI 



CONCLUSIONS AND RECOMMENDATIONS 



As stated previously, optimal control theory states that 
orbital maintenance using a technique in which thrust cancels 
drag, is not optimal. This means that some scheme using 
periodic thrusting must then be optimal. Through the 
complicated techniques of optimal control theory, a thrust 
vectoring scheme is shown to indeed be the optimum. If the 
thrust vector always points along the primer vector, the 
trajectory is optimal. A sub-optimal scheme using fixed-angle 
thrusting and parametric variation is presented here as a 
simplified method of determining the optimality of orbital 
maintenance . 

In each series of tests, minimization of propellent mass 
using fixed-angle thrusting has proven to be inferior to that 
in which thrust is set equal to drag. At first glance these 
results appear to contradict the theory developed by Ross and 
Melton [Ref. 3]. For small perturbations forced-Keplarian 
motion proved to be superior to periodic fixed-angle 
thrusting. As the perturbations increased (indicated by the 
increasing size of the radial band) , it would seem reasonable 
to expect that the difference in fuel consumption between 
these two techniques would increase. However, the results of 
tests described in section A. 4 of Chapter V reveal that for 
large radial bands, the percent difference in propellant mass 
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expended between methods of orbital maintenance using periodic 
thrusting and forced Keplarian motion, actually appears to 
decrease. Based on these results, we may conjecture that the 
percent difference between the two methods tested here will 
continue to decrease until periodic thrusting yields results 
superior to those for forced Keplarian motion. Further testing 
is required before the analytical theory proposed in [Ref. 3] 
may be conclusively verified. 

The problem as presented here is that of fuel-minimization 
during orbital maintenance. Solving the Mayer optimality 
problem derived in Ross [Ref. 3] yields the primer vector. 
This is a very cumbersome method requiring solution of a two- 
point boundary value problem. If periodic thrusting is done 
along the primer vector, fuel will be optimized. This thesis 
has proposed a simpler method using the energy balance of the 
satellite to achieve similar results. Results, however, 
indicate that for small perturbations, f orced-Keplarian motion 
will provide the best results. 

Ultimately, the goal of this thesis is to provide a method 
of fuel-minimization that is practical and may be applied to 
existing systems. Propulsion systems utilizing vectored thrust 
are highly complex and have a correspondingly higher chance of 
failure. Additionally, the extreme size of the perturbations 
required before periodic orbital maintenance would become more 
economical than f orced-Keplarian motion is impractical. 



34 



Based on these conclusions, and with the added knowledge 
that a variable-thrust propulsion system capable of operating 
continuously over the lifetime of a satellite may also be 
impractical, it is recommended that further testing of 
periodic fixed-angle transverse thrusting schemes for small 
perturbations be accomplished. It is recommended that the 
Mayer optimality problem described in Ross [Ref. 3] be solved 
and the primer vector found. The results should then be 
compared to those described in this thesis to determine the 
actual amount of savings acquired through optimization. It is 
possible that the amount of fuel saved may not warrant the 
cost of building a propulsion system capable of periodic 
primer thrusting. Additionally, a comparison of orbital 
maintenance techniques presently in use, with results found in 
this thesis, should be accomplished to determine their 
relative efficiency. 
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APPENDIX A 
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OBJECTIVE: DETERMINATION OF FIXED THRUST ANGLE TO MAINTAIN 
ORBITAL DEVIATION WITHIN A PREDETERMINED 
BANDWIDTH. 



VARIABLE DEFINITIONS: 

X(l) = RADIUS (METERS) 

X (2 ) = RADIAL VELOCITY (METERS PER SECOND) 

X ( 3 ) = THETA (RADIANS) 

X(4) = ANGULAR VELOCITY (RADS PER SECOND) 

XDOT(l) = TIME DERIVATIVE OF X(l) 

XDOT ( 2 ) = TIME DERIVATIVE OF X(2) 

XDOT ( 3 ) = TIME DERIVATIVE OF X(3) 

XDOT (4) = TIME DERIVATIVE OF X(4) 

RO = REFERENCE ORBIT 
D = DRAG (N) 

EO = SPECIFIC ENERGY OF REFERENCE ORBIT 
E = SPECIFIC ENERGY 
MU = GRAVITATIONAL CONSTANT 
M = SPACECRAFT MASS (KG) 

GAMMAR = FLIGHT PATH ANGLE (RADIANS) 
GAMMAD = FLIGHT PATH ANGLE (DEGREES) 

TH = THRUST (N) 

TMAX = BLOWDOWN (MAXIMUM) THRUST 
ALPHAR = THRUST ANGLE (RADIANS) 

CD = COEFFICIENT OF DRAG 

RHOO = REFERENCE ATMOSPHERIC DENSITY 

RHO = CALCULATED ATMOSPHERIC DENSITY 

SPI = SPECIFIC IMPULSE 

V = VELOCITY 

SREF = REFERENCE SURFACE AREA 
GO = GRAVITATIONAL ACCELERATION 
H = INCREMENT OF TIME (STEP SIZE) 

PTI = PRINT TIME INTERVAL (STEP SIZE) 

T = BEGIN TIME 

TF = FINAL (END) TIME 

a = SEMI-MAJOR AXIS 
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C e = ECCENTRICITY 

C 

c 

C START PROGRAM 

PROGRAM ACTORB 
C 

c 

C VARIABLE DECLARATION 

IMPLICIT REAL*8 (A-H, L-Z) 

DIMENSION X ( 4 ) , XDOT ( 4 ) 

c 

c 

C DEFINITION OF CONSTANTS 

PI=3 . 14159265359D+0 
G0=9 . 806D+0 
Tl=0D+0 
N=4 
C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C MAIN PROGRAM 

c 

OPEN (10, FILE=' INIT' , STATUS=' OLD' ) 

OPEN (11, FILE = ' OUT' , STATUS=' NEW' ) 

OPEN (12, FILE=' ORBEL' , STATUS='NEW' ) 

OPEN (13, F ILE=' RAT' , STATUS=' NEW' ) 

C 

READ (10, 1)RO,VO,M,TMAX,T,TF,H,PTI,CD,MU,RHOO,SREF,SPI 
1 FORMAT (13 (/,21X,D13.7) ) 

C 

C 

PRINT*, 'ENTER ALPHA' 

READ*, DEG 

ALPHAR=DEG*P 1/180 . 0D+0 
C 
C 

INDEX=0 

KOUNT=l 

X (1) =R0 

X (2 ) =0D+0 

X (3) =0D+0 

X ( 4 ) =1 . 1 6734 4 9D-3 

E0=M* ( (V0*V0) /2-MU/R0) 

C 
C 

WRITE (11, * 

*ANGLM 
WRITE (11, * 

* I 

WRITE (12, * 



) , ' TIME ORBITS RADIUS VELOCITY ALPHA 

ENERGY TMAX' 

),' (sec) (km) (km/sec) (deg) 

) , ' T a e APOGEE PERIGEE 
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‘PERIOD ' 

WRITE (13, *), ' TIME ORBITS DRAG TH 

* MASS GAMMA' 

C 

C 

C CALCULATIONS 

C 

c 

100 CALL DRAG ( SREF, CD, X, R0 , RHO0 , D, T, V) 

CALL ONOFF ( R0 , X , TMAX , SP I , G 0 , TH , E 0 , GAMMAR, ALPH AR, V, D ) 

CALL DIFFEQ (X, XDOT, MU, D, M, TMAX, ALPHAR, T, TH) 

CALL RK4 (T,X, XDOT, N,H, INDEX, T) 

C 

C 

IF (INDEX .NE. 0) GO TO 100 
C 

c 

C UNIT CONVERSIONS 

C 

R=X ( 1 ) /1000 

SPEED= ( ( ( X ( 2 ) *X(2) ) + (X ( 1 ) *X(4) )**2)**0.5) /1000 

V=SPEED*1000 

M=M- ( (TH*H) / ( SPI *G0 ) ) 

MASS=20000-M 

GAMMAR=ATAN (X (2) / (X ( 1 ) *X (4 ) ) ) 

GAMMAD=GAMMAR* 180. 0D+0/PI 
ANGM= (X ( 1 ) * V*COS (GAMMAR) ) 

ALPHA= ALPHAR* 180 . 0D+0/PI 
ENERGY= ( ( V* V) /2 ) - (MU/X ( 1 ) ) 

ORBITS=T/ 5382 .458 
C 
C 

CALL ORBDAT (ENERGY, ANGM, PI , a, e, APOGEE, PERIGE, PERIOD) 

C 

C 

C OUTPUT 

c 

IF (KOUNT .LT. DNINT (PTI/H) ) GO TO 200 
C 

WRITE (11, 2 ) T, ORBITS, R, SPEED, ALPHA, ANGM, ENERGY, TH 
WRITE (12, 3) T, a, e, APOGEE, PERIGE, PERIOD 
WRITE (13, 4) T, ORBITS, D, TH, MASS , GAMMAD 
C 

2 FORMAT (2X, F7 . 0, IX, F5 . 2 , 3X, F8 . 3 , 2X, F6 . 4, IX, F6 . 1, 2X, 

*F14 . 2, 2X, F12 .2, 2X,F5 .0) 

3 FORMAT (2X, F7 .0, IX, F10 . 3, 2X, F4 . 3, IX, F10 . 3, IX, F10 . 3 , 2x, F10 . 2 ) 

4 FORMAT (2X, F7 . 0, 2X, F5 .2, 2X,F12 . 9, IX, F5 . 0, IX, F15 . 3, IX, F10 . 4) 

C 

C 
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KOUNT=0 

200 KOUNT=KOUNT+l 

IF (TF .GE. T) GO TO 100 
C 

c 

END 

C 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

SUBROUTINE DRAG (SREF, CD, X, R0 , RHO0 , D, T, V) 

IMPLICIT REAL*8 (A-H, L-Z) 

DIMENSION X ( 4 ) , XDOT ( 4 ) 

C 

c 

V= ( (X (2) *X (2) ) + (X (1) *X (4) ) **2) **0 .5 
BETA=2 . 12D-05 

RHO=RHOO*EXP (-BETA* (X (1) -R0) ) 

D = 0 . 5D0 *RHO*CD* SREF*V*V 
C 
C 

RETURN 

END 

C 

C 

ccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccccccccc 
c 

SUBROUTINE DIFFEQ (X, XDOT, MU, D, M, TMAX, ALPHAR, T, TH) 

IMPLICIT REAL*8 (A-H, L-Z) 

DIMENSION X( 4), XDOT (4) 

C 

c 

A=X (1) *X (4) *X (4) 

B=MU/ (X ( 1 ) *X(1) ) 

C= (D/M) * (X (2) / ( ( (X(2) *X(2) ) + (X(1) *X(4) ) **2) **0 .5) ) 

E= (TH/M) * SIN (ALPHAR) 

C 

C 

P=2*X ( 4 ) *X (2 ) /X(l) 

Q=(D/ (X(l) *M) ) * ( (X(l) *X(4) ) / ( ( (X(2) *X(2) ) + (X(l) *X(4) ) **2) **0,5) ) 
S=(TH/ (X(l) *M) ) *COS (ALPHAR) 

C 

C 

XDOT (1) —X (2) 

XDOT (2) =A-B-C+E 
XDOT (3) =X (4) 

XDOT (4)=-P-Q+S 
C 
C 
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9 9 RETURN 

END 

C 

C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

SUBROUTINE RK4 (T, X, XDOT, N, H, INDEX) 

IMPLICIT REAL*8 (A-H, L-Z) 

INTEGER INDEX, I 

DIMENSION X (4) , XDOT (4) , SAVED (4) , SAVEX (4) 

C 

C 

INDEX=INDEX+1 

GO TO (1,2, 3,4), INDEX 

1 DO 10 1=1, N 
SAVEX ( I ) =X ( I ) 

SAVED (I)=XDOT (I) 

10 X ( I ) =SAVEX ( I ) + . 5D0 *H*XDOT ( I ) 

T=T+ . 5D0*H 

RETURN 

C 

C 

2 DO 20 1=1, N 

SAVED (I) =SAVED (I) +2 ,DO*XDOT (I) 

20 X ( I ) =SAVEX ( I ) + . 5D0 *H*XDOT ( I ) 

RETURN 

C 

c 

3 DO 30 1=1, N 

SAVED (I) =SAVED (I) +2 ,D0*XDOT (I) 

30 X (I) =SAVEX (I) +H*XDOT (I) 

T=T+ . 5D0 *H 
RETURN 
C 
C 

4 DO 40 1=1, N 

40 X (I) =SAVEX (I) +H/6 .DO* ( SAVED ( I ) +XDOT ( I ) ) 

INDEX=0 

RETURN 

END 

C 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

SUBROUTINE ORBDAT (ENERGY, ANGM, PI , a, e, APOGEE, PERIGE, PERIOD) 
IMPLICIT REAL*8 (A-H, L-Z) 

C 

C 

e=(l+ (2*ENERGY*ANGM*ANGM/ (3. 9860120813 3D+14* 
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*3. 98 601208 133D+14) ) ) **0.5 
a= (-3 . 986D+14/ (2*ENERGY) ) / 1 0 0 0 
APOGEE=a * (1+e) 

PERIGE=a * (1-e) 

PERIOD=( (2*PI) / (3. 986D+05) **0.5)*(a**1.5) 

C 

C 

RETURN 

END 

C 

C 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

SUBROUTINE ONOFF <R0,X,TMAX, SPI , GO , TH, E0 , GAMMAR, ALPHAR, D) 
IMPLICIT REAL*8 (A-H,L-Z) 

DIMENSION X { 4 ) , XDOT ( 4 ) 

C 

c 

DELTAR=R0-1 0 0 0 

V= ( (X (2) *X(2) ) + (X ( 1 ) *X(4) ) **2) **0 .5 
MU=3 . 986012D+14 
M=2 0 0 0 0 

E=M* ( (V*V) /2- (MU/X ( 1 ) ) ) 

C 

C PRINT*, 'MU=' , MU 

C PRINT*, ' E0 = ' , E0 

C 

IF (TH .EQ. TMAX) GO TO 99 
C 



IF <X(1) .LE. DELTAR) THEN 
C 

IF (E .LT. E0) THEN 
TH=TMAX 

ELSE 

TH=0D+0 

ENDIF 

ELSE 

TH=0D+0 

ENDIF 

C 

c 

99 IF (E .GE. E0) THEN 

TH=0D+0 

ENDIF 

C 

C 

100 RETURN 
END 
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c 



c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
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